home *** CD-ROM | disk | FTP | other *** search
/ SGI Developer Toolbox 6.1 / SGI Developer Toolbox 6.1 - Disc 4.iso / public / pbmplus / pgm / pgmtexture.c < prev    next >
C/C++ Source or Header  |  1994-08-01  |  24KB  |  1,061 lines
  1. /* pgmtxtur.c - calculate textural features on a portable graymap
  2. **
  3. ** Author: James Darrell McCauley
  4. **         Texas Agricultural Experiment Station
  5. **         Department of Agricultural Engineering
  6. **         Texas A&M University
  7. **         College Station, Texas 77843-2117 USA
  8. **
  9. ** Code written partially taken from pgmtofs.c in the PBMPLUS package
  10. ** by Jef Poskanzer.
  11. **
  12. ** Algorithms for calculating features (and some explanatory comments) are
  13. ** taken from:
  14. **
  15. **   Haralick, R.M., K. Shanmugam, and I. Dinstein. 1973. Textural features
  16. **   for image classification.  IEEE Transactions on Systems, Man, and
  17. **   Cybertinetics, SMC-3(6):610-621.
  18. **
  19. ** Copyright (C) 1991 Texas Agricultural Experiment Station, employer for
  20. ** hire of James Darrell McCauley
  21. **
  22. ** Permission to use, copy, modify, and distribute this software and its
  23. ** documentation for any purpose and without fee is hereby granted, provided
  24. ** that the above copyright notice appear in all copies and that both that
  25. ** copyright notice and this permission notice appear in supporting
  26. ** documentation.  This software is provided "as is" without express or
  27. ** implied warranty.
  28. **
  29. ** THE TEXAS AGRICULTURAL EXPERIMENT STATION (TAES) AND THE TEXAS A&M
  30. ** UNIVERSITY SYSTEM (TAMUS) MAKE NO EXPRESS OR IMPLIED WARRANTIES
  31. ** (INCLUDING BY WAY OF EXAMPLE, MERCHANTABILITY) WITH RESPECT TO ANY
  32. ** ITEM, AND SHALL NOT BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL
  33. ** OR CONSEQUENTAL DAMAGES ARISING OUT OF THE POSESSION OR USE OF
  34. ** ANY SUCH ITEM. LICENSEE AND/OR USER AGREES TO INDEMNIFY AND HOLD
  35. ** TAES AND TAMUS HARMLESS FROM ANY CLAIMS ARISING OUT OF THE USE OR
  36. ** POSSESSION OF SUCH ITEMS.
  37. ** 
  38. ** Modification History:
  39. ** 24 Jun 91 - J. Michael Carstensen <jmc@imsor.dth.dk> supplied fix for 
  40. **             correlation function.
  41. */
  42.  
  43. #include "pgm.h"
  44. #include <math.h>
  45.  
  46. #define RADIX 2.0
  47. #define EPSILON 0.000000001
  48. #define BL  "Angle                 "
  49. #define F1  "Angular Second Moment "
  50. #define F2  "Contrast              "
  51. #define F3  "Correlation           "
  52. #define F4  "Variance              "
  53. #define F5  "Inverse Diff Moment   "
  54. #define F6  "Sum Average           "
  55. #define F7  "Sum Variance          "
  56. #define F8  "Sum Entropy           "
  57. #define F9  "Entropy               "
  58. #define F10 "Difference Variance   "
  59. #define F11 "Difference Entropy    "
  60. #define F12 "Meas of Correlation-1 "
  61. #define F13 "Meas of Correlation-2 "
  62. #define F14 "Max Correlation Coeff "
  63.  
  64. #define SIGN(x,y) ((y)<0 ? -fabs(x) : fabs(x))
  65. #define DOT fprintf(stderr,".")
  66. #define SWAP(a,b) {y=(a);(a)=(b);(b)=y;}
  67.  
  68. void results (), hessenberg (), mkbalanced (), reduction (), simplesrt ();
  69. float f1_asm (), f2_contrast (), f3_corr (), f4_var (), f5_idm (),
  70.  f6_savg (), f7_svar (), f8_sentropy (), f9_entropy (), f10_dvar (),
  71.  f11_dentropy (), f12_icorr (), f13_icorr (), f14_maxcorr (), *vector (),
  72.  **matrix ();
  73.  
  74.  
  75. void main (argc, argv)
  76.   int argc;
  77.   char *argv[];
  78. {
  79.   FILE *ifp;
  80.   register gray **grays, *gP;
  81.   int tone[PGM_MAXMAXVAL], R0, R45, R90, angle, d = 1, x, y;
  82.   int argn, rows, cols, bps, padright, row, col;
  83.   int itone, jtone, tones;
  84.   float **P_matrix0, **P_matrix45, **P_matrix90, **P_matrix135;
  85.   float ASM[4], contrast[4], corr[4], var[4], idm[4], savg[4];
  86.   float sentropy[4], svar[4], entropy[4], dvar[4], dentropy[4];
  87.   float icorr[4], maxcorr[4];
  88.   gray maxval, nmaxval;
  89.   char *usage = "[-d <d>] [pgmfile]";
  90.  
  91.     pgm_init( &argc, argv );
  92.  
  93.     argn = 1;
  94.  
  95.     /* Check for flags. */
  96.     if ( argn < argc && argv[argn][0] == '-' )
  97.     {
  98.     if ( argv[argn][1] == 'd' )
  99.         {
  100.         ++argn;
  101.         if ( argn == argc || sscanf( argv[argn], "%d", &d ) != 1 )
  102.         pm_usage( usage );
  103.         }
  104.     else
  105.         pm_usage( usage );
  106.     ++argn;
  107.     }
  108.  
  109.     if ( argn < argc )
  110.     {
  111.     ifp = pm_openr( argv[argn] );
  112.     ++argn;
  113.     }
  114.     else
  115.     ifp = stdin;
  116.  
  117.     if ( argn != argc )
  118.     pm_usage( usage );
  119.  
  120.   grays = pgm_readpgm (ifp, &cols, &rows, &maxval);
  121.   pm_close (ifp);
  122.  
  123.   /* Determine the number of different gray scales (not maxval) */
  124.   for (row = PGM_MAXMAXVAL; row >= 0; --row)
  125.     tone[row] = -1;
  126.   for (row = rows - 1; row >= 0; --row)
  127.     for (col = 0; col < cols; ++col)
  128.       tone[grays[row][col]] = grays[row][col];
  129.   for (row = PGM_MAXMAXVAL, tones = 0; row >= 0; --row)
  130.     if (tone[row] != -1)
  131.       tones++;
  132.   fprintf (stderr, "(Image has %d graylevels.)\n", tones);
  133.  
  134.   /* Collapse array, taking out all zero values */
  135.   for (row = 0, itone = 0; row <= PGM_MAXMAXVAL; row++)
  136.     if (tone[row] != -1)
  137.       tone[itone++] = tone[row];
  138.   /* Now array contains only the gray levels present (in ascending order) */
  139.  
  140.   /* Allocate memory for gray-tone spatial dependence matrix */
  141.   P_matrix0 = matrix (0, tones, 0, tones);
  142.   P_matrix45 = matrix (0, tones, 0, tones);
  143.   P_matrix90 = matrix (0, tones, 0, tones);
  144.   P_matrix135 = matrix (0, tones, 0, tones);
  145.   for (row = 0; row < tones; ++row)
  146.     for (col = 0; col < tones; ++col)
  147.     {
  148.       P_matrix0[row][col] = P_matrix45[row][col] = 0;
  149.       P_matrix90[row][col] = P_matrix135[row][col] = 0;
  150.     }
  151.  
  152.   /* Find gray-tone spatial dependence matrix */
  153.   fprintf (stderr, "(Computing spatial dependence matrix...");
  154.   for (row = 0; row < rows; ++row)
  155.     for (col = 0; col < cols; ++col)
  156.       for (x = 0, angle = 0; angle <= 135; angle += 45)
  157.       {
  158.     while (tone[x] != grays[row][col])
  159.       x++;
  160.     if (angle == 0 && col + d < cols)
  161.     {
  162.       y = 0;
  163.       while (tone[y] != grays[row][col + d])
  164.         y++;
  165.       P_matrix0[x][y]++;
  166.       P_matrix0[y][x]++;
  167.     }
  168.     if (angle == 90 && row + d < rows)
  169.     {
  170.       y = 0;
  171.       while (tone[y] != grays[row + d][col])
  172.         y++;
  173.       P_matrix90[x][y]++;
  174.       P_matrix90[y][x]++;
  175.     }
  176.     if (angle == 45 && row + d < rows && col - d >= 0)
  177.     {
  178.       y = 0;
  179.       while (tone[y] != grays[row + d][col - d])
  180.         y++;
  181.       P_matrix45[x][y]++;
  182.       P_matrix45[y][x]++;
  183.     }
  184.     if (angle == 135 && row + d < rows && col + d < cols)
  185.     {
  186.       y = 0;
  187.       while (tone[y] != grays[row + d][col + d])
  188.         y++;
  189.       P_matrix135[x][y]++;
  190.       P_matrix135[y][x]++;
  191.     }
  192.       }
  193.   /* Gray-tone spatial dependence matrices are complete */
  194.  
  195.   /* Find normalizing constants */
  196.   R0 = 2 * rows * (cols - 1);
  197.   R45 = 2 * (rows - 1) * (cols - 1);
  198.   R90 = 2 * (rows - 1) * cols;
  199.  
  200.   /* Normalize gray-tone spatial dependence matrix */
  201.   for (itone = 0; itone < tones; ++itone)
  202.     for (jtone = 0; jtone < tones; ++jtone)
  203.     {
  204.       P_matrix0[itone][jtone] /= R0;
  205.       P_matrix45[itone][jtone] /= R45;
  206.       P_matrix90[itone][jtone] /= R90;
  207.       P_matrix135[itone][jtone] /= R45;
  208.     }
  209.  
  210.   fprintf (stderr, " done.)\n");
  211.   fprintf (stderr, "(Computing textural features");
  212.   fprintf (stdout, "\n");
  213.   DOT;
  214.   fprintf (stdout,
  215.        "%s         0         45         90        135        Avg\n",
  216.        BL);
  217.  
  218.   ASM[0] = f1_asm (P_matrix0, tones);
  219.   ASM[1] = f1_asm (P_matrix45, tones);
  220.   ASM[2] = f1_asm (P_matrix90, tones);
  221.   ASM[3] = f1_asm (P_matrix135, tones);
  222.   results (F1, ASM);
  223.  
  224.   contrast[0] = f2_contrast (P_matrix0, tones);
  225.   contrast[1] = f2_contrast (P_matrix45, tones);
  226.   contrast[2] = f2_contrast (P_matrix90, tones);
  227.   contrast[3] = f2_contrast (P_matrix135, tones);
  228.   results (F2, contrast);
  229.  
  230.  
  231.   corr[0] = f3_corr (P_matrix0, tones);
  232.   corr[1] = f3_corr (P_matrix45, tones);
  233.   corr[2] = f3_corr (P_matrix90, tones);
  234.   corr[3] = f3_corr (P_matrix135, tones);
  235.   results (F3, corr);
  236.  
  237.   var[0] = f4_var (P_matrix0, tones);
  238.   var[1] = f4_var (P_matrix45, tones);
  239.   var[2] = f4_var (P_matrix90, tones);
  240.   var[3] = f4_var (P_matrix135, tones);
  241.   results (F4, var);
  242.  
  243.  
  244.   idm[0] = f5_idm (P_matrix0, tones);
  245.   idm[1] = f5_idm (P_matrix45, tones);
  246.   idm[2] = f5_idm (P_matrix90, tones);
  247.   idm[3] = f5_idm (P_matrix135, tones);
  248.   results (F5, idm);
  249.  
  250.   savg[0] = f6_savg (P_matrix0, tones);
  251.   savg[1] = f6_savg (P_matrix45, tones);
  252.   savg[2] = f6_savg (P_matrix90, tones);
  253.   savg[3] = f6_savg (P_matrix135, tones);
  254.   results (F6, savg);
  255.  
  256.   sentropy[0] = f8_sentropy (P_matrix0, tones);
  257.   sentropy[1] = f8_sentropy (P_matrix45, tones);
  258.   sentropy[2] = f8_sentropy (P_matrix90, tones);
  259.   sentropy[3] = f8_sentropy (P_matrix135, tones);
  260.   svar[0] = f7_svar (P_matrix0, tones, sentropy[0]);
  261.   svar[1] = f7_svar (P_matrix45, tones, sentropy[1]);
  262.   svar[2] = f7_svar (P_matrix90, tones, sentropy[2]);
  263.   svar[3] = f7_svar (P_matrix135, tones, sentropy[3]);
  264.   results (F7, svar);
  265.   results (F8, sentropy);
  266.  
  267.   entropy[0] = f9_entropy (P_matrix0, tones);
  268.   entropy[1] = f9_entropy (P_matrix45, tones);
  269.   entropy[2] = f9_entropy (P_matrix90, tones);
  270.   entropy[3] = f9_entropy (P_matrix135, tones);
  271.   results (F9, entropy);
  272.  
  273.   dvar[0] = f10_dvar (P_matrix0, tones);
  274.   dvar[1] = f10_dvar (P_matrix45, tones);
  275.   dvar[2] = f10_dvar (P_matrix90, tones);
  276.   dvar[3] = f10_dvar (P_matrix135, tones);
  277.   results (F10, dvar);
  278.  
  279.   dentropy[0] = f11_dentropy (P_matrix0, tones);
  280.   dentropy[1] = f11_dentropy (P_matrix45, tones);
  281.   dentropy[2] = f11_dentropy (P_matrix90, tones);
  282.   dentropy[3] = f11_dentropy (P_matrix135, tones);
  283.   results (F11, dentropy);
  284.  
  285.   icorr[0] = f12_icorr (P_matrix0, tones);
  286.   icorr[1] = f12_icorr (P_matrix45, tones);
  287.   icorr[2] = f12_icorr (P_matrix90, tones);
  288.   icorr[3] = f12_icorr (P_matrix135, tones);
  289.   results (F12, icorr);
  290.  
  291.   icorr[0] = f13_icorr (P_matrix0, tones);
  292.   icorr[1] = f13_icorr (P_matrix45, tones);
  293.   icorr[2] = f13_icorr (P_matrix90, tones);
  294.   icorr[3] = f13_icorr (P_matrix135, tones);
  295.   results (F13, icorr);
  296.  
  297.   maxcorr[0] = f14_maxcorr (P_matrix0, tones);
  298.   maxcorr[1] = f14_maxcorr (P_matrix45, tones);
  299.   maxcorr[2] = f14_maxcorr (P_matrix90, tones);
  300.   maxcorr[3] = f14_maxcorr (P_matrix135, tones);
  301.   results (F14, maxcorr);
  302.  
  303.  
  304.   fprintf (stderr, " done.)\n");
  305.   exit (0);
  306. }
  307.  
  308. float f1_asm (P, Ng)
  309.   float **P;
  310.   int Ng;
  311.  
  312. /* Angular Second Moment */
  313. {
  314.   int i, j;
  315.   float sum = 0;
  316.  
  317.   for (i = 0; i < Ng; ++i)
  318.     for (j = 0; j < Ng; ++j)
  319.       sum += P[i][j] * P[i][j];
  320.  
  321.   return sum;
  322.  
  323.   /*
  324.    * The angular second-moment feature (ASM) f1 is a measure of homogeneity
  325.    * of the image. In a homogeneous image, there are very few dominant
  326.    * gray-tone transitions. Hence the P matrix for such an image will have
  327.    * fewer entries of large magnitude.
  328.    */
  329. }
  330.  
  331.  
  332. float f2_contrast (P, Ng)
  333.   float **P;
  334.   int Ng;
  335.  
  336. /* Contrast */
  337. {
  338.   int i, j, n;
  339.   float sum = 0, bigsum = 0;
  340.  
  341.   for (n = 0; n < Ng; ++n)
  342.   {
  343.     for (i = 0; i < Ng; ++i)
  344.       for (j = 0; j < Ng; ++j)
  345.     if ((i - j) == n || (j - i) == n)
  346.       sum += P[i][j];
  347.     bigsum += n * n * sum;
  348.  
  349.     sum = 0;
  350.   }
  351.   return bigsum;
  352.  
  353.   /*
  354.    * The contrast feature is a difference moment of the P matrix and is a
  355.    * measure of the contrast or the amount of local variations present in an
  356.    * image.
  357.    */
  358. }
  359.  
  360. float f3_corr (P, Ng)
  361.   float **P;
  362.   int Ng;
  363.  
  364. /* Correlation */
  365. {
  366.   int i, j;
  367.   float sum_sqrx = 0, sum_sqry = 0, tmp, *px;
  368.   float meanx =0 , meany = 0 , stddevx, stddevy;
  369.  
  370.   px = vector (0, Ng);
  371.   for (i = 0; i < Ng; ++i)
  372.     px[i] = 0;
  373.  
  374.   /*
  375.    * px[i] is the (i-1)th entry in the marginal probability matrix obtained
  376.    * by summing the rows of p[i][j]
  377.    */
  378.   for (i = 0; i < Ng; ++i)
  379.     for (j = 0; j < Ng; ++j)
  380.       px[i] += P[i][j];
  381.  
  382.  
  383.   /* Now calculate the means and standard deviations of px and py */
  384.   /*- fix supplied by J. Michael Christensen, 21 Jun 1991 */
  385.   /*- further modified by James Darrell McCauley, 16 Aug 1991 
  386.    *     after realizing that meanx=meany and stddevx=stddevy
  387.    */
  388.   for (i = 0; i < Ng; ++i)
  389.   {
  390.     meanx += px[i]*i;
  391.     sum_sqrx += px[i]*i*i;
  392.   }
  393.   meany = meanx;
  394.   sum_sqry = sum_sqrx;
  395.   stddevx = sqrt (sum_sqrx - (meanx * meanx));
  396.   stddevy = stddevx;
  397.  
  398.   /* Finally, the correlation ... */
  399.   for (tmp = 0, i = 0; i < Ng; ++i)
  400.     for (j = 0; j < Ng; ++j)
  401.       tmp += i*j*P[i][j];
  402.  
  403.   return (tmp - meanx * meany) / (stddevx * stddevy);
  404.   /*
  405.    * This correlation feature is a measure of gray-tone linear-dependencies
  406.    * in the image.
  407.    */
  408. }
  409.  
  410.  
  411. float f4_var (P, Ng)
  412.   float **P;
  413.   int Ng;
  414.  
  415. /* Sum of Squares: Variance */
  416. {
  417.   int i, j;
  418.   float mean = 0, var = 0;
  419.  
  420.   /*- Corrected by James Darrell McCauley, 16 Aug 1991
  421.    *  calculates the mean intensity level instead of the mean of
  422.    *  cooccurrence matrix elements 
  423.    */
  424.   for (i = 0; i < Ng; ++i)
  425.     for (j = 0; j < Ng; ++j)
  426.       mean += i * P[i][j];
  427.  
  428.   for (i = 0; i < Ng; ++i)
  429.     for (j = 0; j < Ng; ++j)
  430.       var += (i + 1 - mean) * (i + 1 - mean) * P[i][j];
  431.  
  432.   return var;
  433. }
  434.  
  435. float f5_idm (P, Ng)
  436.   float **P;
  437.   int Ng;
  438.  
  439. /* Inverse Difference Moment */
  440. {
  441.   int i, j;
  442.   float idm = 0;
  443.  
  444.   for (i = 0; i < Ng; ++i)
  445.     for (j = 0; j < Ng; ++j)
  446.       idm += P[i][j] / (1 + (i - j) * (i - j));
  447.  
  448.   return idm;
  449. }
  450.  
  451. float Pxpy[2 * PGM_MAXMAXVAL];
  452.  
  453. float f6_savg (P, Ng)
  454.   float **P;
  455.   int Ng;
  456.  
  457. /* Sum Average */
  458. {
  459.   int i, j;
  460.   extern float Pxpy[2 * PGM_MAXMAXVAL];
  461.   float savg = 0;
  462.  
  463.   for (i = 0; i <= 2 * Ng; ++i)
  464.     Pxpy[i] = 0;
  465.  
  466.   for (i = 0; i < Ng; ++i)
  467.     for (j = 0; j < Ng; ++j)
  468.       Pxpy[i + j + 2] += P[i][j];
  469.   for (i = 2; i <= 2 * Ng; ++i)
  470.     savg += i * Pxpy[i];
  471.  
  472.   return savg;
  473. }
  474.  
  475.  
  476. float f7_svar (P, Ng, S)
  477.   float **P, S;
  478.   int Ng;
  479.  
  480. /* Sum Variance */
  481. {
  482.   int i, j;
  483.   extern float Pxpy[2 * PGM_MAXMAXVAL];
  484.   float var = 0;
  485.  
  486.   for (i = 0; i <= 2 * Ng; ++i)
  487.     Pxpy[i] = 0;
  488.  
  489.   for (i = 0; i < Ng; ++i)
  490.     for (j = 0; j < Ng; ++j)
  491.       Pxpy[i + j + 2] += P[i][j];
  492.  
  493.   for (i = 2; i <= 2 * Ng; ++i)
  494.     var += (i - S) * (i - S) * Pxpy[i];
  495.  
  496.   return var;
  497. }
  498.  
  499. float f8_sentropy (P, Ng)
  500.   float **P;
  501.   int Ng;
  502.  
  503. /* Sum Entropy */
  504. {
  505.   int i, j;
  506.   extern float Pxpy[2 * PGM_MAXMAXVAL];
  507.   float sentropy = 0;
  508.  
  509.   for (i = 0; i <= 2 * Ng; ++i)
  510.     Pxpy[i] = 0;
  511.  
  512.   for (i = 0; i < Ng; ++i)
  513.     for (j = 0; j < Ng; ++j)
  514.       Pxpy[i + j + 2] += P[i][j];
  515.  
  516.   for (i = 2; i <= 2 * Ng; ++i)
  517.     sentropy -= Pxpy[i] * log10 (Pxpy[i] + EPSILON);
  518.  
  519.   return sentropy;
  520. }
  521.  
  522.  
  523. float f9_entropy (P, Ng)
  524.   float **P;
  525.   int Ng;
  526.  
  527. /* Entropy */
  528. {
  529.   int i, j;
  530.   float entropy = 0;
  531.  
  532.   for (i = 0; i < Ng; ++i)
  533.     for (j = 0; j < Ng; ++j)
  534.       entropy += P[i][j] * log10 (P[i][j] + EPSILON);
  535.  
  536.   return -entropy;
  537. }
  538.  
  539.  
  540. float f10_dvar (P, Ng)
  541.   float **P;
  542.   int Ng;
  543.  
  544. /* Difference Variance */
  545. {
  546.   int i, j, tmp;
  547.   extern float Pxpy[2 * PGM_MAXMAXVAL];
  548.   float sum = 0, sum_sqr = 0, var = 0;
  549.  
  550.   for (i = 0; i <= 2 * Ng; ++i)
  551.     Pxpy[i] = 0;
  552.  
  553.   for (i = 0; i < Ng; ++i)
  554.     for (j = 0; j < Ng; ++j)
  555.       Pxpy[abs (i - j)] += P[i][j];
  556.  
  557.   /* Now calculate the variance of Pxpy (Px-y) */
  558.   for (i = 0; i < Ng; ++i)
  559.   {
  560.     sum += Pxpy[i];
  561.     sum_sqr += Pxpy[i] * Pxpy[i];
  562.   }
  563.   tmp = Ng * Ng;
  564.   var = ((tmp * sum_sqr) - (sum * sum)) / (tmp * tmp);
  565.  
  566.   return var;
  567. }
  568.  
  569. float f11_dentropy (P, Ng)
  570.   float **P;
  571.   int Ng;
  572.  
  573. /* Difference Entropy */
  574. {
  575.   int i, j, tmp;
  576.   extern float Pxpy[2 * PGM_MAXMAXVAL];
  577.   float sum = 0, sum_sqr = 0, var = 0;
  578.  
  579.   for (i = 0; i <= 2 * Ng; ++i)
  580.     Pxpy[i] = 0;
  581.  
  582.   for (i = 0; i < Ng; ++i)
  583.     for (j = 0; j < Ng; ++j)
  584.       Pxpy[abs (i - j)] += P[i][j];
  585.  
  586.   for (i = 0; i < Ng; ++i)
  587.     sum += Pxpy[i] * log10 (Pxpy[i] + EPSILON);
  588.  
  589.   return -sum;
  590. }
  591.  
  592. float f12_icorr (P, Ng)
  593.   float **P;
  594.   int Ng;
  595.  
  596. /* Information Measures of Correlation */
  597. {
  598.   int i, j, tmp;
  599.   float *px, *py;
  600.   float hx = 0, hy = 0, hxy = 0, hxy1 = 0, hxy2 = 0;
  601.  
  602.   px = vector (0, Ng);
  603.   py = vector (0, Ng);
  604.  
  605.   /*
  606.    * px[i] is the (i-1)th entry in the marginal probability matrix obtained
  607.    * by summing the rows of p[i][j]
  608.    */
  609.   for (i = 0; i < Ng; ++i)
  610.   {
  611.     for (j = 0; j < Ng; ++j)
  612.     {
  613.       px[i] += P[i][j];
  614.       py[j] += P[i][j];
  615.     }
  616.   }
  617.  
  618.   for (i = 0; i < Ng; ++i)
  619.     for (j = 0; j < Ng; ++j)
  620.     {
  621.       hxy1 -= P[i][j] * log10 (px[i] * py[j] + EPSILON);
  622.       hxy2 -= px[i] * py[j] * log10 (px[i] * py[j] + EPSILON);
  623.       hxy -= P[i][j] * log10 (P[i][j] + EPSILON);
  624.     }
  625.  
  626.   /* Calculate entropies of px and py - is this right? */
  627.   for (i = 0; i < Ng; ++i)
  628.   {
  629.     hx -= px[i] * log10 (px[i] + EPSILON);
  630.     hy -= py[i] * log10 (py[i] + EPSILON);
  631.   }
  632. /*  fprintf(stderr,"hxy1=%f\thxy=%f\thx=%f\thy=%f\n",hxy1,hxy,hx,hy); */
  633.   return ((hxy - hxy1) / (hx > hy ? hx : hy));
  634. }
  635.  
  636. float f13_icorr (P, Ng)
  637.   float **P;
  638.   int Ng;
  639.  
  640. /* Information Measures of Correlation */
  641. {
  642.   int i, j;
  643.   float *px, *py;
  644.   float hx = 0, hy = 0, hxy = 0, hxy1 = 0, hxy2 = 0;
  645.  
  646.   px = vector (0, Ng);
  647.   py = vector (0, Ng);
  648.  
  649.   /*
  650.    * px[i] is the (i-1)th entry in the marginal probability matrix obtained
  651.    * by summing the rows of p[i][j]
  652.    */
  653.   for (i = 0; i < Ng; ++i)
  654.   {
  655.     for (j = 0; j < Ng; ++j)
  656.     {
  657.       px[i] += P[i][j];
  658.       py[j] += P[i][j];
  659.     }
  660.   }
  661.  
  662.   for (i = 0; i < Ng; ++i)
  663.     for (j = 0; j < Ng; ++j)
  664.     {
  665.       hxy1 -= P[i][j] * log10 (px[i] * py[j] + EPSILON);
  666.       hxy2 -= px[i] * py[j] * log10 (px[i] * py[j] + EPSILON);
  667.       hxy -= P[i][j] * log10 (P[i][j] + EPSILON);
  668.     }
  669.  
  670.   /* Calculate entropies of px and py */
  671.   for (i = 0; i < Ng; ++i)
  672.   {
  673.     hx -= px[i] * log10 (px[i] + EPSILON);
  674.     hy -= py[i] * log10 (py[i] + EPSILON);
  675.   }
  676. /*  fprintf(stderr,"hx=%f\thxy2=%f\n",hx,hxy2); */
  677.   return (sqrt (abs (1 - exp (-2.0 * (hxy2 - hxy)))));
  678. }
  679.  
  680. float f14_maxcorr (P, Ng)
  681.   float **P;
  682.   int Ng;
  683.  
  684. /* Returns the Maximal Correlation Coefficient */
  685. {
  686.   int i, j, k;
  687.   float *px, *py, **Q;
  688.   float *x, *iy, tmp;
  689.  
  690.   px = vector (0, Ng);
  691.   py = vector (0, Ng);
  692.   Q = matrix (1, Ng + 1, 1, Ng + 1);
  693.   x = vector (1, Ng);
  694.   iy = vector (1, Ng);
  695.  
  696.   /*
  697.    * px[i] is the (i-1)th entry in the marginal probability matrix obtained
  698.    * by summing the rows of p[i][j]
  699.    */
  700.   for (i = 0; i < Ng; ++i)
  701.   {
  702.     for (j = 0; j < Ng; ++j)
  703.     {
  704.       px[i] += P[i][j];
  705.       py[j] += P[i][j];
  706.     }
  707.   }
  708.  
  709.   /* Find the Q matrix */
  710.   for (i = 0; i < Ng; ++i)
  711.   {
  712.     for (j = 0; j < Ng; ++j)
  713.     {
  714.       Q[i + 1][j + 1] = 0;
  715.       for (k = 0; k < Ng; ++k)
  716.     Q[i + 1][j + 1] += P[i][k] * P[j][k] / px[i] / py[k];
  717.     }
  718.   }
  719.  
  720.   /* Balance the matrix */
  721.   mkbalanced (Q, Ng);
  722.   /* Reduction to Hessenberg Form */
  723.   reduction (Q, Ng);
  724.   /* Finding eigenvalue for nonsymetric matrix using QR algorithm */
  725.   hessenberg (Q, Ng, x, iy);
  726.   /* simplesrt(Ng,x); */
  727.   /* Returns the sqrt of the second largest eigenvalue of Q */
  728.   for (i = 2, tmp = x[1]; i <= Ng; ++i)
  729.     tmp = (tmp > x[i]) ? tmp : x[i];
  730.   return sqrt (x[Ng - 1]);
  731. }
  732.  
  733. float *vector (nl, nh)
  734.   int nl, nh;
  735. {
  736.   float *v;
  737.  
  738.   v = (float *) malloc ((unsigned) (nh - nl + 1) * sizeof (float));
  739.   if (!v)
  740.     fprintf (stderr, "memory allocation failure"), exit (1);
  741.   return v - nl;
  742. }
  743.  
  744.  
  745. float **matrix (nrl, nrh, ncl, nch)
  746.   int nrl, nrh, ncl, nch;
  747.  
  748. /* Allocates a float matrix with range [nrl..nrh][ncl..nch] */
  749. {
  750.   int i;
  751.   float **m;
  752.  
  753.   /* allocate pointers to rows */
  754.   m = (float **) malloc ((unsigned) (nrh - nrl + 1) * sizeof (float *));
  755.   if (!m)
  756.     fprintf (stderr, "memory allocation failure"), exit (1);
  757.   m -= ncl;
  758.  
  759.   /* allocate rows and set pointers to them */
  760.   for (i = nrl; i <= nrh; i++)
  761.   {
  762.     m[i] = (float *) malloc ((unsigned) (nch - ncl + 1) * sizeof (float));
  763.     if (!m[i])
  764.       fprintf (stderr, "memory allocation failure"), exit (2);
  765.     m[i] -= ncl;
  766.   }
  767.   /* return pointer to array of pointers to rows */
  768.   return m;
  769. }
  770.  
  771. void results (c, a)
  772.   char *c;
  773.   float *a;
  774. {
  775.   int i;
  776.  
  777.   DOT;
  778.   fprintf (stdout, "%s", c);
  779.   for (i = 0; i < 4; ++i)
  780.     fprintf (stdout, "% 1.3e ", a[i]);
  781.   fprintf (stdout, "% 1.3e\n", (a[0] + a[1] + a[2] + a[3]) / 4);
  782. }
  783.  
  784. void simplesrt (n, arr)
  785.   int n;
  786.   float arr[];
  787. {
  788.   int i, j;
  789.   float a;
  790.  
  791.   for (j = 2; j <= n; j++)
  792.   {
  793.     a = arr[j];
  794.     i = j - 1;
  795.     while (i > 0 && arr[i] > a)
  796.     {
  797.       arr[i + 1] = arr[i];
  798.       i--;
  799.     }
  800.     arr[i + 1] = a;
  801.   }
  802. }
  803.  
  804. void mkbalanced (a, n)
  805.   float **a;
  806.   int n;
  807. {
  808.   int last, j, i;
  809.   float s, r, g, f, c, sqrdx;
  810.  
  811.   sqrdx = RADIX * RADIX;
  812.   last = 0;
  813.   while (last == 0)
  814.   {
  815.     last = 1;
  816.     for (i = 1; i <= n; i++)
  817.     {
  818.       r = c = 0.0;
  819.       for (j = 1; j <= n; j++)
  820.     if (j != i)
  821.     {
  822.       c += fabs (a[j][i]);
  823.       r += fabs (a[i][j]);
  824.     }
  825.       if (c && r)
  826.       {
  827.     g = r / RADIX;
  828.     f = 1.0;
  829.     s = c + r;
  830.     while (c < g)
  831.     {
  832.       f *= RADIX;
  833.       c *= sqrdx;
  834.     }
  835.     g = r * RADIX;
  836.     while (c > g)
  837.     {
  838.       f /= RADIX;
  839.       c /= sqrdx;
  840.     }
  841.     if ((c + r) / f < 0.95 * s)
  842.     {
  843.       last = 0;
  844.       g = 1.0 / f;
  845.       for (j = 1; j <= n; j++)
  846.         a[i][j] *= g;
  847.       for (j = 1; j <= n; j++)
  848.         a[j][i] *= f;
  849.     }
  850.       }
  851.     }
  852.   }
  853. }
  854.  
  855.  
  856. void reduction (a, n)
  857.   float **a;
  858.   int n;
  859. {
  860.   int m, j, i;
  861.   float y, x;
  862.  
  863.   for (m = 2; m < n; m++)
  864.   {
  865.     x = 0.0;
  866.     i = m;
  867.     for (j = m; j <= n; j++)
  868.     {
  869.       if (fabs (a[j][m - 1]) > fabs (x))
  870.       {
  871.     x = a[j][m - 1];
  872.     i = j;
  873.       }
  874.     }
  875.     if (i != m)
  876.     {
  877.       for (j = m - 1; j <= n; j++)
  878.     SWAP (a[i][j], a[m][j])  
  879.     for (j = 1; j <= n; j++)
  880.       SWAP (a[j][i], a[j][m]) 
  881.       a[j][i] = a[j][i];
  882.     }
  883.     if (x)
  884.     {
  885.       for (i = m + 1; i <= n; i++)
  886.       {
  887.     if (y = a[i][m - 1])
  888.     {
  889.       y /= x;
  890.       a[i][m - 1] = y;
  891.       for (j = m; j <= n; j++)
  892.         a[i][j] -= y * a[m][j];
  893.       for (j = 1; j <= n; j++)
  894.         a[j][m] += y * a[j][i];
  895.     }
  896.       }
  897.     }
  898.   }
  899. }
  900.  
  901. void hessenberg (a, n, wr, wi)
  902.   float **a, wr[], wi[];
  903.   int n;
  904.  
  905. {
  906.   int nn, m, l, k, j, its, i, mmin;
  907.   float z, y, x, w, v, u, t, s, r, q, p, anorm;
  908.  
  909.   anorm = fabs (a[1][1]);
  910.   for (i = 2; i <= n; i++)
  911.     for (j = (i - 1); j <= n; j++)
  912.       anorm += fabs (a[i][j]);
  913.   nn = n;
  914.   t = 0.0;
  915.   while (nn >= 1)
  916.   {
  917.     its = 0;
  918.     do
  919.     {
  920.       for (l = nn; l >= 2; l--)
  921.       {
  922.     s = fabs (a[l - 1][l - 1]) + fabs (a[l][l]);
  923.     if (s == 0.0)
  924.       s = anorm;
  925.     if ((float) (fabs (a[l][l - 1]) + s) == s)
  926.       break;
  927.       }
  928.       x = a[nn][nn];
  929.       if (l == nn)
  930.       {
  931.     wr[nn] = x + t;
  932.     wi[nn--] = 0.0;
  933.       }
  934.       else
  935.       {
  936.     y = a[nn - 1][nn - 1];
  937.     w = a[nn][nn - 1] * a[nn - 1][nn];
  938.     if (l == (nn - 1))
  939.     {
  940.       p = 0.5 * (y - x);
  941.       q = p * p + w;
  942.       z = sqrt (fabs (q));
  943.       x += t;
  944.       if (q >= 0.0)
  945.       {
  946.         z = p + SIGN (z, p); 
  947.         wr[nn - 1] = wr[nn] = x + z;
  948.         if (z)
  949.           wr[nn] = x - w / z;
  950.         wi[nn - 1] = wi[nn] = 0.0;
  951.       }
  952.       else
  953.       {
  954.         wr[nn - 1] = wr[nn] = x + p;
  955.         wi[nn - 1] = -(wi[nn] = z);
  956.       }
  957.       nn -= 2;
  958.     }
  959.     else
  960.     {
  961.       if (its == 30)
  962.         fprintf (stderr, 
  963.                     "Too many iterations to required to find %s\ngiving up\n", 
  964.                      F14), exit (1);
  965.       if (its == 10 || its == 20)
  966.       {
  967.         t += x;
  968.         for (i = 1; i <= nn; i++)
  969.           a[i][i] -= x;
  970.         s = fabs (a[nn][nn - 1]) + fabs (a[nn - 1][nn - 2]);
  971.         y = x = 0.75 * s;
  972.         w = -0.4375 * s * s;
  973.       }
  974.       ++its;
  975.       for (m = (nn - 2); m >= l; m--)
  976.       {
  977.         z = a[m][m];
  978.         r = x - z;
  979.         s = y - z;
  980.         p = (r * s - w) / a[m + 1][m] + a[m][m + 1];
  981.         q = a[m + 1][m + 1] - z - r - s;
  982.         r = a[m + 2][m + 1];
  983.         s = fabs (p) + fabs (q) + fabs (r);
  984.         p /= s;
  985.         q /= s;
  986.         r /= s;
  987.         if (m == l)
  988.           break;
  989.         u = fabs (a[m][m - 1]) * (fabs (q) + fabs (r));
  990.         v = fabs (p) * (fabs (a[m - 1][m - 1]) + fabs (z) + fabs (a[m + 1][m + 1]));
  991.         if ((float) (u + v) == v)
  992.           break;
  993.       }
  994.       for (i = m + 2; i <= nn; i++)
  995.       {
  996.         a[i][i - 2] = 0.0;
  997.         if (i != (m + 2))
  998.           a[i][i - 3] = 0.0;
  999.       }
  1000.       for (k = m; k <= nn - 1; k++)
  1001.       {
  1002.         if (k != m)
  1003.         {
  1004.           p = a[k][k - 1];
  1005.           q = a[k + 1][k - 1];
  1006.           r = 0.0;
  1007.           if (k != (nn - 1))
  1008.         r = a[k + 2][k - 1];
  1009.           if (x = fabs (p) + fabs (q) + fabs (r))
  1010.           {
  1011.         p /= x;
  1012.         q /= x;
  1013.         r /= x;
  1014.           }
  1015.         }
  1016.         if (s = SIGN (sqrt (p * p + q * q + r * r), p)) 
  1017.         {
  1018.           if (k == m)
  1019.           {
  1020.         if (l != m)
  1021.           a[k][k - 1] = -a[k][k - 1];
  1022.           }
  1023.           else
  1024.         a[k][k - 1] = -s * x;
  1025.           p += s;
  1026.           x = p / s;
  1027.           y = q / s;
  1028.           z = r / s;
  1029.           q /= p;
  1030.           r /= p;
  1031.           for (j = k; j <= nn; j++)
  1032.           {
  1033.         p = a[k][j] + q * a[k + 1][j];
  1034.         if (k != (nn - 1))
  1035.         {
  1036.           p += r * a[k + 2][j];
  1037.           a[k + 2][j] -= p * z;
  1038.         }
  1039.         a[k + 1][j] -= p * y;
  1040.         a[k][j] -= p * x;
  1041.           }
  1042.           mmin = nn < k + 3 ? nn : k + 3;
  1043.           for (i = l; i <= mmin; i++)
  1044.           {
  1045.         p = x * a[i][k] + y * a[i][k + 1];
  1046.         if (k != (nn - 1))
  1047.         {
  1048.           p += z * a[i][k + 2];
  1049.           a[i][k + 2] -= p * r;
  1050.         }
  1051.         a[i][k + 1] -= p * q;
  1052.         a[i][k] -= p;
  1053.           }
  1054.         }
  1055.       }
  1056.     }
  1057.       }
  1058.     } while (l < nn - 1);
  1059.   }
  1060. }
  1061.